% 参数初始化
SNR0 = 10; % 初始信噪比（dB），更接近实际场景
N_subcarriers = 128; % 子载波数量
N_symbols = 7; % 符号数量
N_antennas = 2; % 接收天线数量
tau_max = 0.1; % 无线信道的最大延迟扩散（秒）
delta_F = 15e3; % 子载波间隔（Hz）
f_D = 100; % 多普勒频移（Hz）
T_s = 1e-3; % 符号间隔（秒）

% 生成随机信道响应（假设为瑞利衰落）
H_true = (randn(N_subcarriers, N_symbols, N_antennas) + 1j * randn(N_subcarriers, N_symbols, N_antennas)) / sqrt(2);

% 生成DM-RS信号（假设为已知的参考信号）
X_DMRS_OCC = (randn(N_subcarriers, N_symbols) + 1j * randn(N_subcarriers, N_symbols)) / sqrt(2);

% 添加噪声
noise_power = 10^(-SNR0/10); % 噪声功率（线性单位）
noise = sqrt(noise_power/2) * (randn(N_subcarriers, N_symbols, N_antennas) + 1j * randn(N_subcarriers, N_symbols, N_antennas));
Y_RS = H_true .* X_DMRS_OCC + noise; % 接收到的参考信号

% MMSE信道估计
theta = sinc(2 * f_D * T_s * (0:N_symbols-1)) .* sinc(2 * tau_max * delta_F * (0:N_subcarriers-1)');
Phi_inv = 1 ./ (1 + 10^(SNR0/10) * theta);
w = theta .* Phi_inv;
H_est = Y_RS ./ X_DMRS_OCC; % 初次信道估计

% 计算噪声功率
delta = Y_RS - H_est .* X_DMRS_OCC; % 残差
sigma2_hat = mean(abs(delta(:)).^2); % 噪声功率估计

% 计算真实信噪比
X_e = sum(H_est .* X_DMRS_OCC, [1, 2]); % 估计的信号功率
SNR_est = 10 * log10(mean(abs(X_e(:)).^2) / sigma2_hat); % 信噪比估计（dB）

% 第二次信道估计（使用真实信噪比）
Phi_inv_new = 1 ./ (1 + 10^(SNR_est/10) * theta);
H_est_new = Y_RS ./ X_DMRS_OCC; % 更新信道估计

% 计算信道估计的MSE
MSE = mean(abs(H_true(:) - H_est_new(:)).^2); % 均方误差

% 输出结果
fprintf('Estimated noise power: %f\n', sigma2_hat);
fprintf('Estimated SNR: %f dB\n', SNR_est);
fprintf('MSE of channel estimation: %f\n', MSE);

% 绘制时域和频域特性
figure;
subplot(2,1,1);
plot(abs(H_est_new(:,1,1)));
title('Frequency Domain Channel Estimation');
xlabel('Subcarrier Index');
ylabel('Channel Magnitude');

subplot(2,1,2);
plot(abs(fftshift(fft(H_est_new(:,1,1)))));
title('Time Domain Channel Estimation');
xlabel('Time Index');
ylabel('Channel Magnitude');